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ABSTRACT 

The heat capacity of neutron matter is studied over the range of densities and tem¬ 
peratures prevailing in neutron-star crusts, allowing for the transition to a superfluid 
phase at temperatures below some critical temperature T s f and including the tran¬ 
sition to the classical limit. Finite temperature Hartree-Fock-Bogoliubov equations 
(FTHFB) are solved and compared to existing approximate expressions. In particu¬ 
lar, the formula given by Levenfish and Yakovlev is found to reproduce the numerical 
results with a high degree of accuracy for temperatures T ^ T s f. In the non-superfluid 
phase, T ^ T s f, the linear approximation is valid only at temperature T -C Tp n (Tp n 
being the Fermi temperature of the neutron gas) which is rarely the case in the shal¬ 
low layers of the neutron star’s crust. A non-perturbative interpolation between the 
quantal and the classical regimes is proposed here. The heat capacity, conveniently 
parametrized solely in terms of T s f, Tp n , and the neutron number density n n , can be 
easily implemented in neutron-star cooling simulations. 
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1 INTRODUCTION 


Neutron stars (NSs) are formed from the catastrophic gravitational core collapse of stars with a mass M > 8M© (IHaensel et al. 
2007 1. During the first tens of seconds after the collapse, the newly formed proto-neutron star with a radius of about 50 km 
stays very hot with internal temperatures of the order of ~ 10 11 — 10 12 K. But it rapidly cools down and shrinks into an 
ordinary neutron star by emitting neutrinos. As the temperature drops, the outer layer s of the star crystallize into a solid crust, 
whose innermost region is permeated by “free” neutrons (see Charnel & Haensel 20081 . for a review of neutron-star crusts). At 
temperatures T lower than the Fermi temperature Tp n , the neutron liquid becomes degenerate. As the temperature reaches 
some critical temperature T„f, the neutron liquid undergoes a transition to a spin-singlet superfluid phase by forming Cooper 
pairs, like electrons in conventional superconductors. Due to its relatively low neutrino emissivity, the crust of a newly-born 
neutron star cools less rapidly than the core and thus stays hotter. As a result, the surface temperature decreases slowly 
during the first ten to hundred years and then suddenly drops when the cooling wave from the core reaches the surface. 
The e volution of the surface temperature of a young neutron star thus depends essentially on the thermal properties of its 
crust ( Lattimer et al. 19941 : Gnedin et al. 2001 fl. 

While the cooling of very young neutron stars has not been observed yet, the thermal relaxatio n of neutron-star crusts 
has been recently monitored in a few quasi-persistent soft X-ray transients (see Page fe Reddv 20121 . for a recent review). In 
these binary systems, usually called Low-Mass X-Ray Transients (LMXRT), a neutron star accretes matter from a companion 
star during severals years or decades, driving the neutron-star crust out of its thermal equilibrium with the core. After the 
accretion stops, the heated crust relaxes towards equilibrium. While the thermal relaxation of young isolated NSs is mostly 
governed by the thermal properties of the densest layers of the NS crust, the thermal relaxation of LMXRT is expected to be 
very sensitive to the t hermal properties of the shallower layers of the crust, especially near the transition between the inner 
and the outer crusts ( Page fc Reddv 20ld ). In this region of the NS crust, T may be not only higher that T s f but also to 
T Fn so that neutrons can be treated as a classical gas: their heat capacity does not increase linearly with T as in the deeper 
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crustal regions where T s f ^ T -C Tp n , but becomes essentially independent of T. Assuming that neutrons are degenerate in 
all regions of the crust may thus considerably overestimates the crustal heat capacity hence also the thermal relaxation time 


of LMXRT. 


In this paper, we study the neutron contribution to the heat capacity of the neutron-star crust focusing on the region 
close to the ne utron drip trans ition. Although the presence of inhomogeneities may alter the neutron heat capacity (see, e.g., 
Margueron fc Sandulescull2012l l. we shall si mply consider here th at the neutron liquid is homogeneous, as generally assumed in 


neutron-star cooling simulations (see, e.g., Page & Reddy 120121 ). Over the past years, the properties of dilute neutron matter 
have been extensively studied. Microscopic calculations based on realistic nucleon-nucle on potentials an d following different 
many-body approaches tend to yield similar results, at least at low enough densities (Baldo & Burgio 20121 1. On the other 
hand, most of these calculations have been restricted to T = 0, and therefore they cannot be directly applied to neutron-star 
cooling simulations. For this reason, we have employed the finite-temperature Hartree-Fock-Bogoliubov (FTHFB) method, 
which is briefly reviewed in Section [2] The neutron heat capacity is studied in Section [3] for both the classical and degenerate 
regimes. 


2 FINITE-TEMPERATURE HARTREE-FOCK-BOGOLIUBOV METHOD 


The neutron pairing phenomenon is studied here usi ng the self-consistent finite-temperature Hartree-Fock-Bogoliubov ( FTH FB) 
meth od with effective nucleon-nucleon interactions ( Ooodmanlll981 ; Sandulescull2004l : Margueron fc Sandulescull2012l : Pastorel 
2012). This method not only provides a convenient parametrization of microscopic neutron-matter calculations based on re¬ 
alistic potentials at T = 0, but can also be applied to determine in a consistent way the properties of neutron matter at finite 
temperatures including the transition between the quantum and classical regimes . Moreover, the same met hod can also be 
applied to study inhomogeneous nuclear matter like neutron-star crusts (see, e.g., Margueron fc Sandulescul 2012 1. 

A key quantity for characterizing the superfluid phase is the pairing gap function A n (k), k being the wave vector. This 
function, which is related to the binding energy of a pair, depends also on the temperature T (A n vanishing at T ^ T a /) and 
on the neutron number density n n . The pairing gap can be determined by solving the following equations (the units used 
throughout this paper are such that the Boltzmann constant is fee = 1): 
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where £ n {k) denotes the neutron single-particle energies, E n (k) = ( e n {k ) — /i„) 2 + A „(k) 2 the neutron quasi-particle ener¬ 
gies, /in the neutron chemical potential, and v(k, k') the matrix elements of the effective pairing interaction. As in conventional 
superconductors, the effective pairing interaction may not necessarily be the same as the effective interaction which deter¬ 
mines th^single^garticl^gropertiesjn dense matter. For the latter, we use a zero-range effective interaction of the Skyrme 
type ( Skvrmel 19581 : Bender et al. 20031 1 
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where r l: j = n — rj, r = (r t + rj)/2, pij = — ift(V; — Vj)/2 is the relative momentum, &i and frj are Pauli spin matrices, 
and P„ is the two-body spin-exchange operator. Skyrme effec tive interactions we re shown to be well suited for reproducing 
microscopic calculations at both zero and finite temperatures ( Fantina et al. 201 jf . Using Eq. J3|, the neutron single-particle 
energies can be expressed as 
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where Un denotes the (self- consistent) mean potential field felt by neutrons whereas m ^ is the (self-consistent) mean neutron 
effective mass given by I Bender et al. 2003^ 
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m„ being the “bare” neutron mass. Since the FTHFB equations 0-0 depend only on the difference e n {k) — /i„, the pairing 
gap function is actually independent of the potential U n ■ It is therefore convenient to introduce a shifted chemical potential 
v n = Pn — U n . In this way, we have £ n (k) — p n = £n(k) — v n with 
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so that Un disappears from the FTHFB equations. 

As for pairing, we shall also consider zero-range effective interaction of the form 

v(ri,rj) =v n (n n (r))5(r ij ) . (7) 

In this case, the matrix elements of the pairing interaction re duce to v(k , k') = v^lrin). It can thus be seen from Eq. 0 that 
the pairing gap function is independent of k. As shown by Charnel ( 20ld l. the pairing strength v n (n n ) can be directly inferred 
from microscopic calculations of the pairing gap function A n (n„) at T = 0 as follows 
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with 


A(x) = log(16a:) + 2%/l + x — 2 log (l + vT+lr) — 4 . (9) 

The only free parameter is the pairing cutoff ea which needs to be introduced to regularize divergences in Eq. 0 arising 
from the zero-range of the interaction. The pairing functional defined in this way ensures that the solution of the FTHFB 
equations 0-@ coincides exactly with A n at T = 0. With the effective interactions considered here, the gap equations 0 
and © for T <T s f reduce to 
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where E(e) = sjA1 + (e — v n ) 2 , and D(e) denotes the density of neutron single-particle states, defined as 
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3 HEAT CAPACITY OF DILUTE NEUTRON MATTER 

The neutron heat capacity Cy (per unit volume) is defined by ( Landau fe Lifchit3 1980a ) 

8S‘ 


Cy{n n ,T) = T 


dT 


where S n is the entropy density given by 

S"(n n ,T) = - J de [/(e) log/(e) + (1 — /(e)) log(l — /(e))] , 
with /(e) = (1 + exp(E(e)/T)) _1 . Injecting Eq. (11411 into (1131) yields the general expression 
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At high enough density, as in th e core o f a neutron star, neutron pairs are very loosely bound so that A n <C Vn- In the 


weak coupling approximation (Landau & Lifchitz 


, Eq. (11011 is replaced by 
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and Eq. d by u n ~ Tf„, where the Fermi temperature is defined by Tf™ = tn{k = ky n ), see e.g. Eq. (El) , an d where 
/cFn = (37r 2 n„) 1 / 3 denotes the Fermi wave number. Equation (11611 can be further approximated by (Landau & Lifchitz 1980b 1 
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with u = A n (T)/T and the integration has been extended to ±oo. In their seminal paper, Levenfish fc Yakovlev (1994) have 


solved Eq. cm to calculate the neutron heat capacity at densities relevant for neutron star cores, and fitted the numerical 
results with the following analytical expression: 

Cy(n n , T) = Roo(u)Cy^ nor \n„,T), (18) 

where (j™( nor ') j s the heat capacity of non-superfluid neutrons, w hich in the quantum regime of strongly degenerate neutrons, 
T -C Tf-ti, is approximately given by ( Landau fe Lifchitd 1980b 1 
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Figure 1. (Color online) Equation of state (left panel) and effective mass (right panel) in neutron matter for the LNS Skyrme 
functi onal dCao et al.ll 2006l) . The grey area represents the range of equations of state obtained from quantum Monte Carlo calcula¬ 
tions jGari(]ol!r"e7"a!"[|20 1 2l i. See text for details. 


In Eq . (1181) . Rqq is a correction factor introduced to account for the effects of superfluidity and was parametrized as ([Leve nfish fe Yakovlev 
1994h 


Roo{u) = ja 0 + -J (ai) 2 + (a 2 w) 2 j exp (b 0 — yjb\+u- 


( 20 ) 


where u = y/1 — t (co — ci/y/r + ci/t), r = T/T s f. 

Equations ©, ©, © have been widely applied in neutron-star cooling simulations. However, the weak-coupling 
approximation used to obtain these expressions may not be very accurate for the dilute neutron liquid that permeates the 
inner crust of a neutron star, especially in the shallowest region where A„ becomes comparable to u„. Note that in this case, 
the chemical potential v n may be very different from the Fermi energy. It was found, for instance, that the low density regime 
could be located in the BCS-BEC crossover, either close to the u nitary limit, or, i n som e cases, almost entering the BEC 


regime (where the chemical potential becomes negative as T —> 0) (Mar gueron et al. 2007 ). The FTH FB equations (110l) - (lllll 


are however able to describe both the weak and the strong coupling regimes ( Nozieres fc Schmitt-Rinklll985 l. For this reason, 
we have preferred to solve the set of coupled equations © and m without any further approximation instead of Eq. ©. 
The predictions for the heat capacity in the nonsuperfluid and superfluid phases will be discussed separately in Sections 13.II 
and 13.21 respectively. 


3.1 Nonsuperfluid phase: from the classical to the degenerate regimes 


In this section, we first s tudy the heat ca pacity of a non-superfluid neutron liquid. For this purpose, we have considered the 
LNS Skyrme functional ( Cao et al. 200 6). which was fitted to many-body calculations based on the Brueckner method and 


using realistic two- and three-body forces. In particular, the parameters of this functional were adjusted so as to reproduce 
both the equation of state and the effective mass in neutron matter as obtained from many-body calculations. As shown 
in Fig. □ t he equation of state given by this functional is in good agreement with more recent quantum Monte Carlo 
calculations ( Gandolfi et al. 2012! '). 

In Fig[2] (a), the neutron heat capacity calculated from Eq. (1151) with A„ = 0 is compared with the results obtained 
using the linear approximation (1191) . Large deviations can be observed in the low density region where n ra < 10 -3 fm~ 3 or 
equivalently p < 10 12 g cm~ 3 (see, e.g ., Pearson et al.[ (2012l). At high temperatures T O- Tp^, the heat capacity tends to the 


classical limit (Landau & Lifchitz 1980ah 
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Figure 2. (Color online) On the left panel (a): heat capacity of nonsuperfluid neutrons as a function of temperature for different values of 
the density comparing the exact results (symbols) with the linear approximation given by Eq. (119b (dashed lines) and with the improved 
analytical formula given in Eq. J25J (solid lines). On the right panel (b): the correction to the heat capacity induced dy the exchange of 
long-wavelength phonons: solid lines have been calculated using x c = 0.2, while dashed with x c = 0.3. See text for details. 


where 


B 



+ 3v4q a 



(23) 


and Aos = Fo/( 1 + Fq) and Ap s = GW (1 + Go), Fo (resp. Go) being the density density (resp. spin-density spin-density) 
Landau parameter i Bavm fe Pethici3ll991 ) and T c — x c Tf is a cut-off temperature. Eq. (1231) is the lowest order approx i matio n 
considering only the monopolar contribution. Highest order contributions have been calculated by Pethick fe Carneirol ( 1973h . 
The correction A C 7 f' FG ' > normalized to C"J' F G ' > is shown in Fig ,|2l (b) for the three densities (dashed-dotted lines) and for two 

icdflQC 


different choices of the parameter x c (see, e.g.. [Bavin & Pethic k 1 99il l. This correction is shown up to the classical temperature 


T c i since beyond this temperature, the system could not be treated with the same quantal approximation methods, and the 
limit is given by the classical expression <E0>- It is a negligible correction in the case of very diluted neutron matter (less than 
1% in most of the cases). This is mainly due to the absolute value of the Landau Parameters Fo and Go which are small 
compared to 1, for the considered densities. 

In order to assess more precisely the validity of the linear approximation m, we have calculated the relative deviation 


S(T, 


\Cy(T, n n ) — Cy( FG \T,rin)\ 
Cy(T, Tin) 


(24) 


for the range of neutron densities encountered in neutron-star crusts (see, e.g., Pearson et al. 20121 '). Lines of constant 6 in 
a temperature-density diagram are plotted in Fig. [3] Whereas <5 is of the order of a few percents or less at the bottom of 
the neutron-star crust for temperatures T < 10 10 K, errors can be much larger in the shallower layers of the crust near the 
neutron-drip point where n n —> 0. 

The exact results for the neutron heat capacity can be represented by the following simple analytical expression 


G^ (FG - ci) (n„,T)«-n„ 


1 — exp ( — 


T 

Tci 


(25) 


where T c i = 'STf ti /'k 2 is a characteristic temperature delimiting the classical and quantum regimes. This expression yields 
the correct asymptotic behaviours: at low temperatures (T/T c i -C 1) Eg. (1251) reduces to Eg. (1191) . while at high temperatures 
(T/T c i S> 1), Eo. (l25l) tends to the classical limit (HHll . As shown in Fig. [5] Eq. (1251 provides a very accurate interpolation 
between the low and high temperature limits. The classical temperature T c i is also represented in Fig. [3] showing that at 
T ss T c i , the relative deviation 5(T, n n ) amounts to ~ 30%. 

Note that in Eq. ([25]). the effect of the neutron-neutron interactions are embedded in the neutron effective mass appearing 
in the expression of the Fermi temperature Tp n . However, rW ~ m n at densities below n n ~ 10~ 2 fm -3 therefore the 
expression (1351) is almost independent of the choice of the functional over this range of densities. 

While in this paper, we are mostly interested in neutron matter, Eq. (|25l) can be easily adapted to evaluate the heat 
capacity of any other fermion species, e.g. electrons. 
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n Q [fm 3 ] 

Figure 3. (Color online) Temperature as a function of density corresponding to fixed values of the error function S(T,n n ). The dashed 
line represents the characteristic temperature T c i appearing in Eq. <25t and delimiting the classical from the quantum regime. See text 
for details. 



Figu re 4. (Color online) 1 So neutron pairing gap as a function of density as obtained from Eq. (126ft (solid line ) and from dCao et all 
l20()d ) (symbols). The gray band is extracted from taking the most recent ab — initio results of Ref. ifCezerlis fc Carlsonlboioh . 


3.2 Superfluid phase 


Let us now study the heat capacity in the superfiuid regime and the transition to the normal phase at T = T s f. In our 
calculations, we us e an effective pairing interac tion which is adjusted to mimic more microscopic intera ctions based on 
Brueckner method ( Brueckner fc Levinson 19551 ). We have adopted the ones obtained by Cao et al. ( 20061 ) which includes 
medium polarization effects beyond the mean-field. The reason for calibrating our effective pairing interaction to these results 
is that they were determined using the same Brueckner approach as the one we have used to fix the mean-field and the 
effective mass, see Section I3T1 In this way, the particle-ho l e and the particle-particle channels are consistently defined. 

The mi croscopic pairing gaps calculated bv lcao et al. ( 20061 ) at T = 0 can be conveniently parametrized by the following 
expression ( Kaminker et alj 2001) 


A n(T = 0) = Q(kmax ~ fc.F„)Ao 


k Fn (fcl-'n ^2) 

k 2 Fn + fc l (k Fn - k 2 ) 2 + fcf 


(26) 

fci, & 2 , &3 are given in Tab.Q] 


in which Q(x) denotes the Heaviside unit step function. The values of the p ara meter s Ao, k„ 

As shown in Fig. |4j expression (1261) can well reproduce the gap by Cao et al. ( 20061 )._ 

In the literature several ab-initio calculations of the pairing gap have been performed ( Gezerlis fe Carlson 20ld) . They 











































Heat capacity of low density neutron matter: from quantum to classical regimes 7 


Table 1. Values of the parameters entering into Eq. lt26ll . 


k-max [fm 1] 

A 0 [MeV] 

fcl [fm 3 ] 

&2 [fm x ] 

[fm J ] 

1.37 

3.37968 

0.556092 

1.38236 

0.327517 



T/T 

C 


Figure 5. (Color online) Evolution of the 1 So neutron pairing gap as a function of temperature for some representative value of the 
density of the system expressed in [fm~ 3 ]. The dots represent the exact FTHFB results, while the solid line is obtained using Eq. 12711 . 


are included in the gray band in Fig© There are some discrepancies among them, as indicated by the thickness of the band, 
but the parametrization of the pairing gap that we consider (solid green line) remains within this band. 

Injecting the gap expression (1261) into the definition of the effective pairing force ©, we have solved the FTHFB Eqs. ©- 
© and we have determined the evolution of the pairing gaps as a function of temperature. The value of the pairing cutoff 
can be arbitrarily chosen. As shown in Fig© the gap decreases with increasing temperature, and eventually vanishes at the 
critical te mperature T s f . The temperature dependence of the pairing gap can be accurately represented by the following 
expression (Goriely 19961 ) 


l (T) = A n (T = 0)\jl-( — 


0(T — T sf ). 


(27) 


for T ^ T 3 f. Th e critica l temper ature can be obtained from the pairing gap at zero temperature from the well-known 
relation (Landau & Lifchitz 1980bl ) 


Tsf = 


exp(£) 


A„(T = 0), 


(28) 


where £ ss 0.577 is the Euler-Mascheroni constant. Equations (12611 . (12711 and (1281) provide a convenient parametrization of the 
evolution the pairing gaps as a function of both temperature and density. 

From the solutions of the FTHFB Eqs. ©-©, we have calculate d the heat capacity u sing Eq. (1151) and compared the 
results with those obtained using the analytical formula proposed by L evenfish fc Yakovlev (119941 ). As shown in Fig. [6] this 
formula reproduces very w ell the FTHF B resu lts in all cases and for the domain of temperature ranging from 0 to few MeV. 

The formula of Levenfish fc Yakovlev ( 1994 ) can be further improved by readjusting the numerical parameters of the Rqo 
form factor. As shown in Fig© the new formula for the heat capacity 


Cv(n n ,T) « -n n R 00 (T) 


1 — exp I — 


T 

Tfi 


(29) 


provides a very food fit to the exact FTHFB results both in the quantal and classical regimes. 


4 CONCLUSIONS 

The heat capacity of dilute neutron matter has been calculated in the framework of the finite-temperature Hartree-Fock- 
Bogoliubov method using an energy density functional fitted to microscopic calculations at zero temperature. Although the 
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Figure 6. (Color online) Neutron heat capacity: FTHFB results (dots) are compared to the original expression ILevenfish &; Yakovlevl 
|l99J), see Eg. <181 1. (solid line) and to the new expression <2511 with refitted parameters given in Tab[2] (dashed line). The values of the 
density are expressed in [fm —3 ]. 



LY 

present 

ao 

0.4186 

-0.7294 

ai 

1.007 

2.331 

d2 

0.501 

0.806 

bo 

1.456 

0.743 

h 

1.456 

0.743 

7 

2.5 

2.5 

Co 

1.456 

-2.689 

Cl 

0.157 

-3.800 

C2 

1.764 

0.826 


Table 2. Values of the parameters appearing in Eg. <201 . LY refers to the original values given bv lLevenfish fe Yakovlevl d 19941 1. 


heat capacity is found to be well reproduced by the empirical expressions of L evenfish fe Yakovlev (|1994 ) at high enough 
density, it differs substantially in the dilute limit corresponding to the shallow region of the inner crust of a neutron star. In 
particular, the use of Eg. (11811 overestimates the heat capacity of nonsuperfluid neutrons, hence also the thermal relaxation 
time of neutron- star crusts. Indeed, n onsuperfluid neutrons provide the main contribution to the heat capacity of the neutron- 
star inner crust ([Gnedin et al.,1120011 ). We have proposed a new analytical expression obtained by combining Eqs. (11811 . (1201) . 
([25]) . that tends to the correct limits both at T = 0 and at T —> +oo. This expression also takes into account the superfluid 
transition at T = T s f, the critical temperature T s / being given by Eqs. (ESI) and ([28]. 

T he presence of neutron-proton clusters in the inner crust of a neutron star may impact the heat capacity of the neutron 


i ne presence of neutron-proton clusters m tnc inner crust ot a neutron star may impact tnc neat capacity ol tne ncutroi 
liquid ( Pizzochero et al. 2002 ; SandulesciJ 2004 ; Monrozeau et al. 2007 ; Fortin et al. 2Q10l : Charnel et al. 201(1 : Pastord 20121 ) 


an effe ct which we haven’t considered in this study. While this approximation was shown to be very accurate in the normal 
phase ( Charnel et al. 20091 ). it may induce larger errors in the superfluid phase. In most regions of the inner crust, the pairing 
gap hence also the critical temperature are reduced as compared to pure neutron matter calculation s, bu t thi s effect can be 


easily included in the present expression for Cy by simply renormalizing T sf (see, e.g., Ichamel et aklboiQ ). As a matter 


of fact, this reduction of T s f lies within the theoretical uncertainties of microscopic calculations in pure neutron matter. On 
the other hand, spatial inhomoge nei ties c an drasti cally change the phase diagram of the neutron liquid near the neutron- 
drip transition ( Margueron fe Khan 20121 : Pastord 2012). The effect of non-uniform matt er could be parameterized as an 
additional correction to the uniform matter expressions, as performed by Fortin et al. 1 20ld ). The simple expression proposed 
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in this work, valid in all regimes from classical’s to superfluid’s ones, could efficiently be used to determine accurately the 
aforementioned corrections. Further studies are therefore needed to determine the neutron heat capacity in these layers of the 
star, but the present work provides an improved reference for more refined calculations taking into account the presence of 
nuclear clusters in the crust. 
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